### Alizade, Dancygier, Ditlmann
### "National Penalties Reversed"
### Replication Code 
### Figure A4
### For questions, contact jalizade@princeton.edu

# setup
rm(list = ls())
setwd("C:/Users/Jey/Dropbox/WZB/NaturalizationExperiment/Submission/JOP/replication_JOP/data")
library(foreign)
library(lmtest)
library(sandwich)
library(ggplot2)
dat <- read.dta("data_experimental.dta")

# convert treatment variable from second experiment to factor
dat$e2_treat[dat$e2_treat=="NA"] <- NA
dat$e2_treat <- factor(dat$e2_treat, levels=c("Vote Intention (Canadian)", "Control (Turkish)", "Vote Intention (Turkish)", "Integration Problems (Turkish)"))

### Figure A4a ###

# create data frame with frequencies by outcome
freq_a <- c(table(dat$e1_legal_dummy[dat$e1_response==1]), 
            table(dat$e1_help_dummy[dat$e1_response==1]), 
            table(dat$e1_posaff_dummy[dat$e1_response==1]))
x_a <- as.numeric(names(freq_a))
out_a <- rep(c("Information", "Help", "Positive Affect"), each=2)
means_a <- rep(c(mean(dat$e1_legal_dummy[dat$e1_response==1 & !is.na(dat$e1_response)]),
                 mean(dat$e1_help_dummy[dat$e1_response==1 & !is.na(dat$e1_response)]),
                 mean(dat$e1_posaff_dummy[dat$e1_response==1 & !is.na(dat$e1_response)])), 
                 each=2)
df_a <- cbind.data.frame(freq_a, x_a, means_a, out_a)
df_a$out_a <- factor(df_a$out_a, levels=c("Information", "Help", "Positive Affect"))

# plot
plot_a <- ggplot(data=df_a, aes(x=x_a, y=freq_a))
plot_a <- plot_a + geom_bar(stat="identity")
plot_a <- plot_a + facet_wrap(~ out_a, scales = "free_x")
plot_a <- plot_a + geom_vline(aes(xintercept = means_a), linetype="dashed", color="red")
plot_a <- plot_a + scale_x_continuous(breaks=c(0,1)) 
plot_a <- plot_a + xlab("At least one category present") + ylab("Frequency")
plot_a


### Figure A4b ###

# create data frame with frequencies by outcome
freq_b <- c(table(dat$e2_legal_dummy[dat$e2_response==1]), 
            table(dat$e2_help_dummy[dat$e2_response==1]), 
            table(dat$e2_posaff_dummy[dat$e2_response==1]))
x_b <- as.numeric(names(freq_b))
out_b <- rep(c("Information", "Help", "Positive Affect"), each=2)
means_b <- rep(c(mean(dat$e2_legal_dummy[dat$e2_response==1 & !is.na(dat$e2_response)]),
                 mean(dat$e2_help_dummy[dat$e2_response==1 & !is.na(dat$e2_response)]),
                 mean(dat$e2_posaff_dummy[dat$e2_response==1 & !is.na(dat$e2_response)])), 
               each=2)
df_b <- cbind.data.frame(freq_b, x_b, means_b, out_b)
df_b$out_b <- factor(df_b$out_b, levels=c("Information", "Help", "Positive Affect"))

# plot
plot_b <- ggplot(data=df_b, aes(x=x_b, y=freq_b))
plot_b <- plot_b + geom_bar(stat="identity")
plot_b <- plot_b + facet_wrap(~ out_b, scales = "free_x")
plot_b <- plot_b + geom_vline(aes(xintercept = means_b), linetype="dashed", color="red")
plot_b <- plot_b + scale_x_continuous(breaks=c(0,1)) 
plot_b <- plot_b + xlab("At least one category present") + ylab("Frequency")
plot_b


### Figure A4c ###

# create data frame with frequencies by outcome
freq_c <- c(table(dat$e1_legal_sum[dat$e1_response==1]), 
            table(dat$e1_help_sum[dat$e1_response==1]), 
            table(dat$e1_posaff_sum[dat$e1_response==1]))
x_c <- as.numeric(names(freq_c))
out_c <- c(rep("Information", 4), rep("Help", 4), rep("Positive Affect", 3))
means_c <- c(rep(mean(dat$e1_legal_sum[dat$e1_response==1 & !is.na(dat$e1_response)]), 4),
             rep(mean(dat$e1_help_sum[dat$e1_response==1 & !is.na(dat$e1_response)]), 4),
             rep(mean(dat$e1_posaff_sum[dat$e1_response==1 & !is.na(dat$e1_response)]), 3))
df_c <- cbind.data.frame(freq_c, x_c, means_c, out_c)
df_c$out_c <- factor(df_c$out_c, levels=c("Information", "Help", "Positive Affect"))

# plot
plot_c <- ggplot(data=df_c, aes(x=x_c, y=freq_c))
plot_c <- plot_c + geom_bar(stat="identity")
plot_c <- plot_c + facet_wrap(~ out_c, scales = "free_x")
plot_c <- plot_c + geom_vline(aes(xintercept = means_c), linetype="dashed", color="red")
plot_c <- plot_c + xlab("Number of relevant coding categories present") + ylab("Frequency")
plot_c


### Figure A4d ###

# create data frame with frequencies by outcome
freq_d <- c(table(dat$e2_legal_sum[dat$e2_response==1]), 
            table(dat$e2_help_sum[dat$e2_response==1]), 
            table(dat$e2_posaff_sum[dat$e2_response==1]))
x_d <- as.numeric(names(freq_d))
out_d <- c(rep("Information", 4), rep("Help", 4), rep("Positive Affect", 3))
means_d <- c(rep(mean(dat$e2_legal_sum[dat$e2_response==1 & !is.na(dat$e2_response)]), 4),
             rep(mean(dat$e2_help_sum[dat$e2_response==1 & !is.na(dat$e2_response)]), 4),
             rep(mean(dat$e2_posaff_sum[dat$e2_response==1 & !is.na(dat$e2_response)]), 3))
df_d <- cbind.data.frame(freq_d, x_d, means_d, out_d)
df_d$out_d <- factor(df_d$out_d, levels=c("Information", "Help", "Positive Affect"))

# plot
plot_d <- ggplot(data=df_d, aes(x=x_d, y=freq_d))
plot_d <- plot_d + geom_bar(stat="identity")
plot_d <- plot_d + facet_wrap(~ out_d, scales = "free_x")
plot_d <- plot_d + geom_vline(aes(xintercept = means_d), linetype="dashed", color="red")
plot_d <- plot_d + xlab("Number of relevant coding categories present") + ylab("Frequency")
plot_d


